Orbital currents in extended Hubbard models of high-T c cuprates 
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Motivated by the recent report of broken time-reversal symmetry and zero momentum magnetic scattering 
in underdoped cuprates, we investigate under which circumstances orbital currents circulating inside a unit cell 
might be stabilized in extended Hubbard models that explicitly include oxygen orbitals. Using Gutz wilier pro- 
jected variational wave functions that treat on an equal footing all instabilities, we show that orbital currents 
indeed develop on finite clusters, and that they are stabilized in the thermodynamic limit if additional interac- 
tions, e.g. strong hybridization with apical oxygens, are included in the model. 



Despite intensive efforts in the last twenty years, the 
physics of high Tc superconductors remains largely mysteri- 
ous [1]. This is especially true of the pseudogap phase of un- 
derdoped cuprates, for which various explanations have been 
put forward ranging from preformed superconducting pairs 
[5] to the existence of orbital currents with [6] or without [7,8] 
broken translational symmetry. The latter case has received 
considerable attention due on one hand to recent neutron ex- 
periments [9] indicating the presence of magnetic moments, 
compatible with the translation invariant pattern of currents 
predicted by Varma [7], and on the other hand to Kerr effect 
measurements [10] showing evidence of time-reversal sym- 
metry breaking. 

Current (flux) phases have been first proposed for the 
single-band Hubbard model [2, 5, 11, 12], but have been 
found unstable by slave bosons and numerical calculations. 
Interestingly, orbital currents were also found to be relevant in 
two-band system [3? ]. In ladder models, where the existence 
of such phases can be checked in a more controlled way, it 
was found that somewhat special interactions, more complex 
than local ones, are needed to stabilize them [13, 14]. The re- 
sulting phases break the translational symmetry of the lattice, 
leading to a staggered flux pattern. Similar staggered patterns 
were advocated as a potential explanation of the pseudogap 
phase [6] (DDW phase). 

To stabilize flux phases that do not break the translational 
symmetry, it seems worthwhile to go beyond the single band 
Hubbard model and to consider the so-called three-band Hub- 
bard model in which oxygen orbitals are explicitly taken into 
account. While this model has been tested for supercon- 
ducting instabilities early on, pointing to some differences 
with the single band model [15], the investigation of its flux 
phase instabilities has only started quite recently. The one- 
dimensional (ladder) case has first been investigated at half 
filling, when the system is essentially an insulator [16]. More 
recently, a low energy analysis for the three band model on a 
ladder has been performed, and showed that in a certain range 
of doping flux phases were indeed stabilized [17]. These 
phases exhibit a 2k f staggered order parameter, quite natu- 
ral in one-dimension, for the currents, but with a symmetry 
different from that of the DDW. In two dimensions, a mean- 
field analysis [7] of the three band model has suggested the 



existence of translationally invariant current patterns when the 
Cu-0 nearest-neighbor repulsion is strong enough. This result 
has not been confirmed by exact diagonalizations [18, 19] on 
small clusters of an effective t—J model. However the map- 
ping of the three-band Hubbard model on this t—J model can 
only be justified in the limit of very large oxygen on-site en- 
ergy, a situation not realized in the cuprates [20]. Moreover, 
because of the three atoms per unit cell, the exact diagonal- 
izations have been performed for clusters with only a few unit 
cells, and the relevance of the results for the thermodynamic 
limit is far from guaranteed, in particular because the filling 
that was considered on this small cluster x = 12.5% is lead- 
ing to a polarized ground- state with finite momentum, which 
is not representative of the physics on large scales. There- 
fore, further investigations of the three-band Hubbard model 
are clearly called for. 

In this paper we perform a variational Monte Carlo (VMC) 
investigation of the three band Hubbard model based on a 
Gutzwiller projected wave function that allows for the pos- 
sibility of orbital currents. This provides an unbiased method, 
free from numerical limitations even for large system sizes, 
for which current instabilities are treated on an equal foot- 
ing with other instabilities. We find that on intermediate sys- 
tem sizes, a current flux phase circulating between copper and 
oxygens is stabilized. This phase has the same symmetry (62, 
see Fig. 1) as the phase found in the mean field solution [7]. 
Other symmetries or phases that break the translational sym- 
metry are much higher in energy. However, as the system 
size gets larger, the energy gain decreases strongly, making 
it unclear whether such a phase would survive in the thermo- 
dynamic limit. We propose modification of the Hamiltonian 
that takes into account apical oxygens or three-body terms and 
which strongly stabilizes such current patterns. 

The three band model is defined by the Hamiltonian: 
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where tij stand for the hopping matrix elements between 
nearest-neighbor Cu-0 pairs ((i, j)), of magnitude td p , and 
between next-nearest neighbor 0-0 pairs of mag- 

nitude tpp, while A p , Ud, U p and Vd p denote the atomic en- 
ergy of the 0-p orbitals, the on-site repulsions in the Cu-d 
and 0-p orbitals, and the nearest-neighbor repulsion between 
Cu-d and 0-p orbitals. A realistic set of parameters found by 
LDA calculations [20, 21, 22] is U d = 10.5 eV and U p = 4 
eV, \t dp \ = 1.3 eV and \t pp \ = 0.65 eV, A p = 3.5 eV, and 
Vd p = 1.2 eV, and, unless specified otherwise, these values 
are used in the following. Note that we work in hole nota- 
tions. The phase factor that comes from the hybridization of 
thep—d orbitals [34] gives to the a non homogeneous sign. 

Around each copper atom, there is one Cu-O-Cu plaquette 
with three minus signs. Moreover, the product of the hopping 
signs around all the Cu-O-Cu plaquette is —1 in hole nota- 
tions. Interestingly, a simple gauge transformation involving 
a double copper unit-cell leads to = — \Uj | [35]. 

This model has already been investigated with variational 
Monte Carlo [23, 24], but the wave functions used did not al- 
low for current instabilities. By contrast, the wave function we 
consider throughout this work is constructed from the ground 
state of the Hofstadter-like mean-field hamiltonian: 
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where Xij > 0,%, A v p ar and hi are real variational parame- 
ters. In this expression, the sum runs over nearest and next- 
nearest neighbor pairs, and c stands for p or d orbitals de- 
pending on the site. The local magnetic field allows to 
consider antiferromagnetism. The 16 variables Xij an d % 
are independent within one copper unit-cell. To deal with 
staggered order we multiply a subset of the variational pa- 
rameters by -1. When 6^ ^ 0, the order parameters are 
associated with an external flux which leads to the circula- 
tion of the holes. In this case the Green function is complex: 
(c\cj) = \(c\cj)\ exp(i^j) We treat the correlations with a 
spin and charge Jastrow factor: 
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where all and v^_j\ are considered as free variational 

parameters. 

We are mainly interested in the charge current observable. 
The conservation of the density: 




Srii he r- n -r 

C (hj) 



St 



(4) 



leads to the current operator [36] on a link: = 

Uj I i c l c j ) I sm (<fiij )• However the current conservation law 

cr 

is not satisfied by our variational wave function since it is not 
an exact eigenstate of (1). The mean-field current, defined 



as J™ := 

a 

quantity, but in general J and Jm f need not be oriented in 
the same way. In particular we find that the 0-0 currents have 
different signs for J and Jm f • The phase e zeij thus gives the 
direction of J MF , but not the direction of J in general, which 
has to be computed explicitly. 

In order to overcome this difficulty and to impose the cur- 
rent conservation at each vertex of the lattice, we apply on the 
wavefunction an additional complex Jastrow factor [25]: 
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The kinetic energy of the extended wavefunction is (T) = 

tij I ( c L c jcr) I cos (Oij + a j ~ a i)- We emphasize that 

the variables are in general not able to cancel the flux 6ij . 
The minimization of the variational parameters introduced in 
(2) and in (3) is performed using a stochastic minimization 
procedure [26, 27] in which the parameters of the uncorrelated 
part of the wavefunction and the Jastrow parameters are mini- 
mized at the same time. This method allows one to deal with a 
large number of parameters since the gradients are calculated 
all at the same time during a simulation. The new parameters 
are then calculated using the obtained gradients, and the pro- 
cedure is iterated until convergence. At each step, the param- 
eters {ai} of (5) are determined by finding the ground state 
of a classical 2D XY spin model. Once our wave function is 
optimized, we measure the physical observables, normalized 
by the number of copper atoms in the lattice. To benchmark 
our wave function, we have compared its properties with those 
obtained by exact diagonalizations for 10 holes on a small 8 
copper lattice with periodic boundary conditions (25% hole 
doping), with very encouraging results (see Table I). Details 
will be presented elsewhere [28]. 

We now turn to larger lattices, using open boundary condi- 
tions to avoid the frustration of the ai variables. We have per- 
formed calculations on lattices ranging from 16 to 64 copper 
sites. The gain in energy for the best antiferromagnetic (SDW) 
instability and for flux phases are compared in Fig. 1 [37] Note 
that these orders are somewhat exclusive: considering both si- 
multaneously does not lead to any detectable gain in energy. 
At zero doping we find that the Neel magnetic long-range or- 
der is stabilized. By introducing doping, we considered as a 
first approximation only the Q = (ir, 7r) pitch vector for the 
spin density wave, even though the pitch vector is most likely 
doping dependent [30] or other instabilities like stripes can oc- 
cur [23]. Nevertheless, the long-range correlations contained 
in the Jastrow factor allow for a decent treatment of the spin 
correlations. Indeed for our best variational w ave functi on, 
the magnetic order parameter M = lim r ^ 00 yj (SfSf +r )) is 

w 60% of the classical value. Using this wavefunction as 
a guiding function for the fixed node calculations, we find 
a slightly higher magnetic order (62%). These values com- 
pare well with the 60% obtained by quantum Monte Carlo in 
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wave function 


Etot 


Tdp 


Tpp 


u d 


Up 


A p 


V d p 


variance 


Lanczos 


-1.13821 


-3.10036 


-0.79666 


0.26737 


0.08398 


1.77545 


0.63201 





Jastrow wf 


-1.0775(1) 


-3.06068 


-0.83073 


0.26176 


0.08197 


1.83466 


0.64076 


0.018 


1 Lanczos step 


-1.1153(1) 


-3.14070 


-0.83715 


0.26559 


0.08736 


1.86495 


0.64469 


0.018 


Fixe node/ Jastrow 


-1.1112(5) 






0.26966 


0.08708 


1.77314 


0.63177 


0.001 



TABLE I: Variational energies of the different wavefunctions compared with exact diagonalizations (Lanczos) on an 8-copper lattice with 10 
holes and S z = 0. We show the total energy (Etot), the kinetic energy of the copper-oxygen links (Td P ) and of the oxygen-oxygen links (T pp ), 
the on-site repulsion energy of the d (Ud) and p (U p ) orbitals, the expectation value of the charge gap operator (A p ), and the expectation value 
of the Coulomb repulsion between the d and p orbitals (Vdp). The results obtained by applying one Lanczos step (the best variational ansatz) 
and the fixed node approximation [29] on a simple Jastrow wavefunction are also shown (the calculation of the off-diagonal operators Tdp 
and Tpp within the fixe-node approximation is more involved and is beyond the scope of the present study). The comparison shows that our 
wavefunction captures quite well the low energy physics of the model. 
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FIG. 1 : Energy of the best orbital current wavefunction (top) and of 
the best spin density wave wavefunction (bottom). The reference was 
taken as the projected Fermi sea (full Jastrow projector). Inset: the 
symmetry of the best current phase found (62 pattern in the notations 
of [7]). Although the orbital current phase is stable on an intermedi- 
ate size 16 copper lattice, the energy gain is strongly reduced when 
the size increases. To the accuracy of the calculation, we did not 
find any energy optimization on a 100 copper lattice for the orbital 
current instability. 



the one-band Heisenberg model [31]. Note that, in the three- 
band Hubbard model, the magnetic instability is strongly de- 
pendent on the oxygen-oxygen hopping integral t pp . Despite 
the good agreement for the order parameter, our variational 
treatment overestimates the magnetic instability, which only 
disappears for 13% doping, instead of the experimentally ob- 
served x = 2% for the cuprates. VMC tends indeed to over- 
estimate the stability of SDW phases because the alternating 
magnetization allows to avoid double occupancy in the uncor- 



rected part of the wavefunction. The presence of magnetic 
order clearly costs kinetic energy, but it does a better job than 
a pure local Gutzwiller projection, which reduces the kinetic 
energy much more dramatically. 

Let us now turn to the current instability. On intermedi- 
ate system sizes such as the 16 copper lattice, we find that 
an orbital current phase with a symmetry 62 (see the inset of 
Fig. 1) is stabilized. Our wave function thus shows a tendency 
to flux phase that does not exist in the corresponding one band 
model. However the gain in energy strongly decreases as the 
size increases and seems, within the accuracy of our calcu- 
lation, to vanish in the thermodynamic limit. Taken literally, 
this suggests that the orbital current phase is not stable, and 
that the system only has short range correlations. However we 
of course cannot rule out that some fine tuning of the param- 
eters could stabilize this phase, or that the energy gain would 
be much smaller than our statistical error. Two points should 
be emphasized. First, regardless of the size of the system, 
we find consistently that the 62 symmetry is the one with the 
lowest variational energy. Other symmetries such as the 61 
phase or the DDW patterns are unstable within our variational 
approach in the whole range of doping. Second, varying the 
Cu-0 interaction seems to have no major effect on the sta- 
bility of the current patterns, in contrast with what one could 
have expected on the basis of the mean-field solution [7]. Our 
study demonstrates that for an instability towards long-range 
current correlations to develop, the relative sign of the transfer 
integrals around a loop is a crucial parameter. Indeed, when 
the sign of t pp is reversed, J and Jm f are oriented in the same 
direction, and the current pattern is stable. In the mean field 
solution, this change of sign is induced by the correction to 
the bare t p d coming from the decoupling of the Cu-0 repul- 
sion term. 

Given our variational results for the three-band Hubbard 
model of Cu02 planes, and the experimental evidence, it is 
natural to check if other terms not included into this simple 
version of a multi-band model could produce such an effect 
and stabilize orbital currents. One obvious possibility is to in- 
clude correlated hopping, as already done successfully for the 
single-band Hubbard model. This indeed leads to a strong in- 
crease of the tendency to develop orbital currents [28]. Note 
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FIG. 2: Circulation of current around a p x —p z —p y plaquette (trian- 
gles) and around a p x — d^ Z 2- r 2 —p y plaquette (squares) measured in 
our best variational Ansatz for the eight band Hubbard model. The 
phenomenological parameter 77 renormalizes the amplitudes of the 
out-of-plane transfer integrals. The symmetry of the orbital current 
pattern is 62 like: there are two out-of-plane current loops in the 
upper pyramid (a) and two current loops in the copper plane (b). 
Finally, the current pattern in the lower pyramid (not shown) is ob- 
tained by a horizontal mirror reflection of the upper pyramid. The 
calculations were done on a 32 Cu lattice at x — 0.125 hole doping. 



however that the magnitude of these terms is difficult to as- 
sess. A second and interesting possibility is to generalize the 
model by including apical oxygens. Indeed, the relative signs 
of the hopping around purely oxygen plaquettes allow a pri- 
ori for orbital currents. We have thus repeated the calculation 
for Cu-0 layers including apical oxygens above and below 
each Cu atom as well as the Cu-d 3z 2_ r 2 orbitals. The bare pa- 
rameters have been taken from LDA [32] performed with the 
crystal structure of the insulating parent compound, and the 
analysis has been carried out as a function of the distance be- 
tween the copper and the apical oxygen, modelled by a renor- 
malization parameter 77 that multiplies the hopping integrals 
involving the apical oxygens. As expected from the signs of 
the hopping, a calculation on a 32-Cu lattice (192 sites) with 
x = 0.125 hole doping confirms that orbital currents involv- 
ing the apical oxygen or the Cu-d 3z 2_ r 2 indeed develop ac- 
cording to the pattern of Fig. 2. These currents are quite small 
for the bare values of the parameters, but they strongly in- 
crease to reach significant values above 77 = 1.5. Interestingly 
enough, the current circulating in the p x — p z — p y plaquette 
leads to a tilted moment, which would provide a natural ex- 
planation for the out of plane moment that was observed in 
the neutron experiments [9] . Whether the structural changes 



induced by doping on the position of the apical oxygens re- 
ported by some authors can be large enough to produce such a 
strong renormalization of the hopping integrals remains to be 
seen [33]. 
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